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Abstract: 

We describe a simple analytic representation of the ArH 2 potential energy sur- 
face which well reproduces the results of extensive ab initio electronic structure 
calculations. The analytic representation smoothly interpolates between the disso- 
ciated H 2 and strong bonding limits. In the fitting process, emphasis is made on 
accurately reproducing regions of the potential expected to be important for high 
temperature (ca. 3000 K) collision processes. Overall, the anisotropy and H 2 bond 
length dependence of the analytic representation well reproduce the input data. 
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I. Introduction 


In this paper we present a new potential energy surface (PES) for ArH 2 - This 
PES is specifically tailored for use in high temperature (ca. 3000 K) energy transfer 
studies involving dissociation and recombination of H 2 . This potential will be used 
to help understand the recombination process 

H + H + M-+H 2 +M. 

The recombination rates for M = H 2 [1], Ar [2] and H 2 O [2] have been computed 
as a function of temperature using potential energy surfaces from the literature 
[3-5]. The result is that within a factor of 1.4 over the temperature range 1000 - 
5000 K, the recombination rate is independent of the third body M. This contrasts 
with the prevailing experimental interpretation [6] which assigns the efficiencies 
5:1:0. 3 for H 2 O, H 2 and Ar at 2000 K. There are several potential sources for this 
disagreement: inaccuracies in the potential energy surfaces used, limitations in the 
dynamics/kinetics calculations (see Ref.[l] for a description of the techniques we 
use to compute the recombination rate given a PES), or experimental uncertainties. 
Certainly of the three third bodies, the PES for M=Ar was the least well known. 
In this paper we give a PES of comparable quality to the potential energy surfaces 
for the other third bodies we have studied. 

When constructing a PES, many different strategies are possible. One can com- 
bine data from various sources, both theoretical and experimental, and attempt to 
produce a function which reproduces this information. This procedure is hampered 
by the possibilities of inconsistency between the various input data. An alternate 
approach, which is used here, is to just consider data from one source, namely high 
quality ab initio electronic structure calculations. This strategy can be criticized 
in that certain features of the PES, such as the weak van der Waals minimum, are 
obtainable more easily from experimental analysis than from theoretical calcula- 
tions. However, in our opinion, this drawback is more than counterbalanced by the 
consistency of the input data which should give rise to a PES which may more accu- 
rately reflect the shape of the accurate PES over wide regions of configuration space. 
We also anticipate that the features of the PES likely to be given most accurately 
by the electronic structure calculations, namely regions of moderate repulsion, will 
be predominately responsible for governing the outcome of the high temperature 
energy transfer processes which are of primary interest to us. 
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In fitting a PES for which a wide range of different interactions are possible, 
one has to choose between a simple function which can be expected to be accurate 
only in the mean or a more complicated function with many parameters. For a 
simple function, it is often easy to ensure that it is well behaved in extrapolated 
regions, while for more complicated functions it may be very difficult to ensure that 
extrapolation to geometries not included in the fitting procedure will give rise to 
physically reasonable results. In the present work we begin by using a very simple 
function to represent the Ar + H 2 interaction and then carefully modify it to ensure 
that it is well behaved globally yet accurately reproduces the input data. This is 
achieved using only eight parameters optimized to fit 71 energy points. 

II. Electronic structure calculations. 

The Ar basis set starts with a (17sl2p6d4f) primitive set. The two outermost s, 
p, and d and the outermost f primitives are left uncontracted and the inner functions 
axe contracted using an atomic natural orbital (ANO) scheme [7] to [5s4p2dlf]. In 
addition a (2s2pldlf) set of even tempered diffuse functions is added. The resulting 
basis is denoted as [5+1+ls 4+1+lp 2+1+ld 1+lf] + (2s2pldlf). 

The H basis set starts with a (8s6p4d) primitive set. The outermost s and p 
primitives axe left uncontracted and the inner functions are contracted based on 
ANO’s to [3s2pld]. This basis is denoted as [3-fls 2+lp Id]. 

The energy calculations were designed to give an accurate description of the 
process H + H + Ar —* H 2 + Ar. The most important nondynamical electron 
correlation effects were taken into account by means of complete active space, self 
consistent field (CASSCF) calculations, and the dynamical correlation effects were 
estimated using the approximately size-extensive averaged coupled pair functional 
(ACPF)[8] method. The CASSCF configurations correspond to a single configura- 
tion description of Ar and a two-electron two active orbited description of H 2 . The 
full CASSCF reference space was used in the ACPF calculations. All calculations 
were carried out using the MOLECULE-SWEDEN [9] program system on the Ames 
Research Center Advanced Computation Facility CRAY Y-MP/864. 

The ANO basis for H 2 was optimized at the equilibrium bond length ( R e ) [7]. 
This makes the basis set superposition error much smaller at R e than at other H 2 
bond lengths. To help correct for this systematic error, we corrected all our energies 
by means of the Boys and Bernardi function counterpoise method[9.1]. Specifically, 
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we computed the following energies: 

P ArH2 ((?), the energy of the ArH 2 supermolecule at the ACPF level at geom- 
etry G , 

P H2 (G), the energy of the H 2 molecule at the all single and double excitation 
configuration interaction (SDCI) level at geometry G including the Ar one- 
electron basis functions, 

E Al (G), the energy of the Ar atom at the ACPF level at geometry G including 
the H atom one-electron basis functions. 

Let G°° be a geometry approximating Ar not interacting with H 2 but having the 
same H 2 bond length as in geometry G. In practice, not interacting means any two 
atom-atom distances greater than 30 clq. Then we take the energy of interaction 
between Ar and H 2 to be 

int 

It should be noted that, since the SDCI and ACPF methods are not exactly size- 
extensive, it is important to use Eq.(l) rather than E Ai ^ 2 (G) — E H2 (G) — E At (G ) 
to compute the interaction energy since E AlH2 (G°°) ^ E Ai (G °° ) + E H2 (G°°). 

In Tables I and II we give the values of P mt computed in the course of this 
work. The geometries are specified by the Jacobi coordinates r, 12, and %, where r 
is the distance between the Ar and the center of mass of the H 2 , R is the H 2 bond 
length, and % is the Ar-center of mass H 2 - H angle. 

III. Fitting the PES 

We will optimize the parameters (JP) in our fitting function [y tnl ((7, P)] by 
minimizing the least squares function 

HP) = ZO-E^Gi) - F i "‘(G i ,.P)]/® i } ! . (2) 

i 

An important aspect of the present work is the choice of the weights Wi. Pre- 
vious work [10] has shown strong correlations between the force along the diatom 
bond evaluated at the classical turning point for relative translation and vibrational 
energy transfer rates, thus we wish to concentrate our interest to regions where 
V int ~ E re i- At 3000 K, k&T = 10 -2 E h. Also u >h for H 2 is about 2 x 10 -2 Ph, 
so we want to focus in most closely on energies in the range 5 — 50 mPh* We will 


(G) = E AlH2 {G)-E AlJl2 (G 00 )+E n2 (G°°)-E n2 (G)+E Al {G co )-E AT {G). (1) 
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identify the weights as an acceptable error in the fit for each point. In the most 
interesting region, a 5% error is reasonable, thus we take 



This is a parabola in % error as a function of the order of magnitude of the energies, 
centered on lOmE'h, except at very weakly interacting geometries, where we let the 
acceptable error be as large as O.limEh* With this weighting scheme, reporting the 
root mean square deviation between the fit and the input data is not meaningful 
— the important quantity is which is unitless. A value equal to one means 

that on the average, we have attained our goal in fitting the points. 

The first step in our fitting process is the determination of an ArH poten- 
tial curve, V ArH . We will use this in a zeroth order representation of the ArH 2 
interaction energy, i.e. 


yint _ V ArH( R ^ + V AiH( R _^ ( 4 ) 

where R± are the two ArH distances, i.e . 

R ± = [r 2 + (R/2) 2 ± rR cos x ] 3 • ( 5 ) 

This simple form will give physically reasonable results for all geometries provided 
yAiK - g we jj behaved. We choose to parameterize the potential curve as 

frArH ffArH /^ArH 

F ArH (X) = 4exp(— bX] - — J 6 ~ + j 4 )2 - (x 2 + <f 2 ) 5 ’ ^ 

where C ArH — C^, rH are taken from Ref. [11]. The three parameters A, b, and d 
axe determined by nonlinear least squares using the weights of Eq.(3) and the data 
given in Table I. The final value of <j> for this fit is 0.23 with maximum weighted 
difference of 0.45. The final values of the fitting parameters are given in Table III 
where they are labeled Ao , &o and do • With these parameters, F ArH is well behaved 
for all atom-atom distances. 

We now turn to geometries where the two H atoms are interacting with the Ar. 
We staxt with the simple sum of pairwise interactions, Eq.(4). This simple form 
has several deficiencies. We correct for these deficiencies by introducing additional 
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dependence on the coordinates r, R : and % other than that implicit in Eq.(5). In 
introducing the additional coordinate dependence, one has considerable flexibility. 
However, since in the limit R goes to oo, Eq.(4) is exact, we require that the 
corrections to Eq.(4) vanish in this limit. The easiest way to enforce this limit is to 
have the corrections depend only on R and not on r or x- We thus considered an 
intermediate step in the fitting process, namely fitting the ab initio data for fixed 
bond lengths. This will yield parameters fully optimized as a function of R . In the 
final step, we will fit the dependence of the parameters on R. The intermediate 
step involves fewer points and parameters than the final step and provides excellent 
initial guesses for nonlinear parameters. It will also tell us whether or not it is 
possible to obtain a satisfactory fit with only R dependent parameters. In the 
present case, we will see that it will be possible to obtain accurate fits with this 
restriction. 

The first correction to Eq.(4) concerns the long range attractive term. The 
leading contribution from /{R^. + d 6 ) expanded in terms of the Jacobi co- 

ordinates is 2 C^ rH /r 6 + 0(r“ 8 ) which can be compared with the more accu- 
rate form C'^' rH2 (f2)[l + Tg 1112 (R)P 2 (cosx)]/r 6 + 0(r -8 ). Our hypothesis requires 
that we neglect the P 2 term, however this is not a bad approximation since 
r£ rH2 (E) ~ 10% for R near R e [13]. We thus scale the long range part of Eq.(4) 
by the factor g(R) = (7^ rH2 (f2)/2(7^ rH and we fit the data of Varandas [13] to 
g(R) = {/3 exp[— a(R — J?o)] + l} -1 with a = 2 ajj” 1 , (3 = 0.437 and Rq — 1.449ao. 
As R becomes large, g(R) becomes unity and we recover the ArH long range inter- 
action. 

The next deficiency we address is the anisotropy of the sum of pairwise inter- 
actions, e.g. the % dependence. Using Eq.(4) with the scaled long range potential 
results in too much anisotropy. We interpret this as a manifestation of the migra- 
tion of charge density toward the center of the H 2 bond, thus we introduce the shift 
/?, i.e. we replace R± in Eq.(5) with 


R± = 


7,2 + ( ^2 ~ ) ±r (' R ~^ > ) cos ^] 


( 7 ) 
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thus we take 


V int {r,R, X ,A,b,d,p) = A[ex V (-bR + )+exp{-bR-)] 

- g(R){C$ lH [(#[. + d 6 )' 1 + (R1 + d 6 )" 1 ] + C 8 ArH [(R\ + d 4 )" 2 + (R 4 _ + d 4 ) -2 ] 
+ C'lV 11 [(R-l + d 2 )- 5 + (R- + d 2 )" 5 ] }• 

( 8 ) 

Finally we optimize the four parameters A , b, d , and p for each value of R. The value 
of 4> for each value of R is given in Table IV under the column marked optimum. A 
satisfactory fit is obtained for all H 2 bond lengths. 

To complete our fit, we must parameterize A(i2), b(R ), d(R ), and p(R). We 
first consider p(R)- Since this parameter is based upon the idea of charge flow, 
we will tie it to a molecular property which reflects this, namely the quadrupole 
moment. Specifically we take 

p(R) = R{ 1 - [1 + q(R)/8R 2 ]-*}, (9) 

where q(R) is the quadrapole moment of H 2 . We will use the fit of Ref. [3] to 
generate q(R) at all R. This contains the single parameter 5. Originally, we tried to 
parameterize p(R) as the displacement of a fixed point charge required to reproduce 
q(R), but we were not successful. Equation (9) represents the displacement required 
to reproduce q(R) for a R dependent point charge with R dependence leading to a 
constant discriminant in the quadratic equation for which p{R) is a root. 

The remaining functions are represented as an expansion in terms of even 
tempered exponentials, i.e. 

^max , X 

X(R)= X i [-exp(axR)] i , (10) 

t=0 

where X = A, b, or d, with the i = 0 coefficient fixed by the fit to the ArH 
interaction. In our final fit we use i msX)J 4 = 2 and im»x,6 = fmax,<2 = 1- The 
parameters for the fit are given in Table III and the weighted errors as a function of 
R in Table IV. The maximum weighted error is 2.7 and the number of geometries 
with weighted errors greater than one is five. With these parameters, the fit seems 
to be globally well behaved. 

IV. Discussion 
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We first consider the ArH interaction. The fit is quite good, however it does 
not well reproduce the empirical ArH potential from Ref.[ll]. See Fig. 1. The 
present ArH potential curve has van der Waals minimum parameters iZmjn = 7.13 ao 
Vjnin — *-0.102 mE^ compared to Rmin — 6.68 a 0 Fmin = —0.175 mi?h for the Tang- 
Toennies potential including up to C is- The Tang-Toennies potential is also much 
more repulsive at distances less than i?min* The reason for this discrepancy is not 
clear, for improving the theoretical calculations by both expanding the one-electron 
basis set and changing the electron correlation treatment does not significantly 
change the present result [12]. 

We can also compare to empirically determined values of the functions in the 
expansion 

V i7lt (r,R = 1.449 a 0 ,%) = v 0 (r) + P 2 (cos %)u 2 (r) (11) 

given in Ref. [14]. This is done in Fig. 2 where we also give v 0 and v 2 computed 
from the ab initio data and the current fit. These functions for the fit do not change 
significantly between R = 1.401 and R = 1.449, so it is meaningful to compare 
the functions for the two distances. The comparison between the v 0 term is very 
similar to the ArH potential: the fit well reproduces the ab initio data and the 
minimum for the Tang-Toennies function is deeper and closer in: iZmin = 6.77 ao 
and Vjnin = —0.232 mE^ versus iZmin = 7.10 a 0 and FmWi = —0.157 m£h for the 
fit. In contrast, for the v 2 term, the Tang-Toennies function fits the ab initio data 
better than the fit in the well region. The failure of the fit to reproduce the ab 
initio result is presumably due to the fact that the leading long range anisotropic 
contribution is proportional to r -8 . However, since v 2 is not large, the absolute 
magnitudes of the deviations are not large. 

Now consider the consequence of stretching the H 2 bond. The most interesting 
orientation to consider is that for % = 7r/2, the T-shaped geometries. In Fig. 3 we 
compare the fit to the ab initio data for repulsive geometries. On the whole we see 
very good agreement between the fit and the data. This particular orientation is 
interesting because the ab initio data and the fit predict that at iZ c , V tni will increase 
as R increases whereas a simple pairwise additive model would give the reverse 
trend. It is the inclusion of the R dependent shift p(R) which is primarily responsible 
for the correct sign of the slope for the present fit. Eventually as iZ increases, the 
slope changes sign and the fit appears to do a reasonable job at predicting the 
position of the sign change. However, in this figure, larger deviations appear for 
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R = 3 than for the other geometries. This is also shown by the weighted errors in 
Table IV. This is probably an indication that a more sophisticated treatment of the 
anisotropy is required. 

In the course of the present work we considered several variants of the present 
fit. All of these were attempts to improve the description of the anisotropy of the 
PES. An obvious place to start is to replace the function g{R ) used to scale the 
dispersion terms with the function of Ref.[13] which includes x dependence and a 
more accurate R dependence. We also considered introducing P 2 ( cos x) anisotropy 
into the functions A(R), b(R), d(R) and 6(R). By doing so, it was possible to 
substantially reduce the weighted errors for R = 1 and 1.401, but the improvement 
for R = 3 and 5 was much more modest. These fits also yielded more complicated 
R dependence for A(R) etc and so it was very hard to produce functions which gave 
rise to a globally well behaved PES. On the whole, the present fit was outstanding 
in its ability to accurately predict the global shape of the PES. 

IV. Conclusions 

In this paper we have presented a simple yet accurate function to interpo- 
late/extrapolate the interaction energy of Ar+H 2 . In conjunction with an accurate 
H 2 potential curve, such as that from Ref. [3], this defines a global PES suitable 
for studying energy transfer and dissociative collisions. The function does not well 
reproduce empirical estimates of the van der Waals well region, but does reproduce 
reasonably well the ab initio data in that region. The force along the H 2 bond is 
well represented globally without recourse to switching functions. 
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Table I. ArH interaction energies. 


^ArH(«o) 

E int (m E h ) 

F ArH (m E h ) 

2. 

215.740 

205.625 

3. 

50.853 

49.295 

4. 

11.278 

11.505 

5. 

1.912 

1.885 

6. 

0.119 

0.110 

7. 

-0.095 

-0.101 

8. 

-0.072 

-0.078 

9. 

-0.040 

-0.044 

10. 

-0.021 

-0.024 
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Table II. ArH2 interaction energies. 




jjjint 

( m £ h ) 


K°o) 

R = 1 a 0 

R = 1.401 

CO 

II 

05 

R = 5 

3 . 

158.061 

x = o 
182.810 

535.157 

a 

4 . 

25.707 

32.480 

73.251 

522.077 

5 . 

3.202 

4.546 

16.575 

90.337 

6 . 

0.089 

0.219 

2.757 

22.709 

8 . 

- 0.094 

- 0.129 

- 0.218 

0.533 

10 . 

- 0.025 

- 0.035 

- 0.082 

- 0.103 

3 . 

144.256 

X = tt /4 
154.744 

151.256 

150.373 

4 . 

23.772 

28.272 

41.471 

62.092 

5 . 

2.981 

3.937 

8.873 

18.126 

6 . 

0.098 

0.190 

1.175 

4.016 

8 . 

- 0.082 

- 0.112 

- 0.185 

- 0.088 

10 . 

- 0.020 

- 0.031 

- 0.062 

- 0.068 

3 . 

133.271 

X = v / 2 
136.783 

97.310 

29.719 

4 . 

22.015 

24.734 

23.799 

8.183 

5 . 

2.765 

3.404 

4.233 

1.444 

6 . 

0.097 

0.164 

0.340 

- 0.016 

8 . 

- 0.077 

- 0.097 

- 0.145 

- 0.119 

10 . 

- 0.021 

- 0.027 

- 0.045 

- 0.037 


a Very repulsive so not computed. 
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Table III. ArH 2 PES parameters (in atomic units). 


X 

X 0 


x 2 

a x 

A 

18.06 

58.63 

-90.61 

0.7662 

b 

1.653 

0.3874 


0.7662 

d 

2.339 







5 = 0.2369 
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Table IV. Root mean square errors for ArH 2 PES fit. 


R 

Optimum 

Fit 

1.0 

0.20 

0.32 

1.401 

0.29 

0.46 

3.0 

0.85 

0.98 

5.0 

0.55 

0.82 

all 

0.54 

0.70 
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Figure captions: 

Figure 1: The ArH potential curve from ab initio calculations (symbols), from the 
present fit , Eq. (6) (solid line), and from Tang and Toennies [11]. Note the switch 
between a linear and log scale at 0.1 mE^. 

Figure 2: The first two terms in a Legendre expansion of the ArH 2 PES. Circles, 
t>o from ab initio points; squares V 2 from ab initio points; solid line fit; dashed fine, 
the fit of Tang and Toennies [14]. Note the switch between a linear and log scale at 
0.1 mEh- 

Figure 3: The dependence of the potential energy on the H 2 bond length for x = 
7r/2. The symbols are the ab initio points with squares for r = 3 ao> circles for 
t — 4, triangles for r = 5, and crosses for r = 6. The solid lines are the fit at the 
corresponding values of r. 
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Ar+H2 

chi=pi/2 


□ r- 3.000 





